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Abstract: We provide estimates of glacier mass changes in the High Mountain Asia (HMA) area from April 
2002 to August 2016 by employing a new version of gravity solutions of the Gravity Recovery and Climate 
Experiment (GRACE) twin-satellite mission. We found a total mass loss trend of the HMA glaciers at a 
rate of —22.17 (+1.96) Gt/a. The largest mass loss rates of —7.02 (+0.94) and —6.73 (£0.78) Gt/a are found 
for the glaciers in Nyaingentanglha Mountains and Eastern Himalayas, respectively. Although most glaciers 
in the HMA area show a mass loss, we find a small glacier mass gain of 1.19 (£0.55) and 0.77 (£0.37) Gt/a 
in Karakoram Mountains and West Kunlun Mountains, respectively. There is also a nearly zero mass balance 
in Pamirs. Our estimates of glacier mass change trends confirm previous results from the analysis of 
altimetry data of the ICESat (Ice, Cloud and Land Elevation) and ASTER (Advanced Spaceborne Thermal 
Emission and Reflection Radiometer) DEM (Digital Elevation Model) satellites in most of the selected 
glacier areas. However, they largely differ to previous GRACE-based studies which we attribute to our 
different post-processing techniques of the newer GRACE data. In addition, we explicitly show regional 
mass change features for both the interannual glacier mass changes and the 14-a averaged seasonal glacier 
mass changes. These changes can be explained in parts by total net precipitation (net snowfall and net 
rainfall) and net snowfall, but mostly by total net radiation energy when compared to data from the ERA5- 
Land meteorological reanalysis. Moreover, nearly all the non-trend interannual mass changes and most 
seasonal mass changes can be explained by the total net radiation energy data. The mass loss trends could 
be partly related to a heat effect due to increased net rainfall in Tianshan Mountains, Qilian Mountains, 
Nyaingentanglha Mountains and Eastern Himalayas. Our new results for the glacier mass change in this 
study could help improve the understanding of glacier variation in the HMA area and contribute to the 
study of global change. They could also serve the utilization of water resources there and in neighboring 
areas. 
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1 Introduction 


In High Mountain Asia (HMA), the glaciers are mostly distributed in Tianshan Mountains 
(Tianshan), Pamirs, Hindu Kush Mountains (Hindu Kush), Karakoram Mountains (Karakoram), 
Himalayas and Nyainqentanglha Mountains (Nyainqentanglha). The total glacier area in HMA 
accounts for about 23% of the global mountain glacier area (Gardner et al., 2013) (Fig. 1). 
Particularly, glacier mass balances (Cogley et al., 2011), also called glacier mass changes in the 
space-gravimetry community (Kääb et al., 2012), are very sensitive to long-term global warming. 
Understanding glacier mass changes is critical for the utilization of water resources in HMA and 
adjacent areas and the study of global change (Brun et al., 2017; Pritchard et al., 2019). For this 
purpose, different survey techniques have been developed and used. 
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Fig. 1 Trend rates of monthly gravity changes (in equivalent water height, EWH) in the High Mountain Asia 
(HMA) and adjacent areas from April 2002 to August 2016, derived from the Gravity Recovery and Climate 
Experiment (GRACE) release-6 solution from the Institute of Geodesy of Technische Universität Graz, Austria 
(ITSG). Mascon 1 covers the glaciers in Nyaingentanglha Mountains (Nyainqentanglha), mascons 2 to 3 in Eastern 
Himalayas, mascons 4 to 5 in Western Himalayas, mascon 6 in Karakoram Mountains (Karakoram), mascon 7 in 
Hindu Kush Mountains (Hindu Kush), mascon 8 in Pamirs, mascons 9 to 12 in Tianshan Mountains (Tianshan), 
mascon 13 in West Kunlun Mountains (West Kunlun), and mascon 14 in Qilian Mountains (Qilian), respectively. 


Glaciological methods can be used to evaluate glacier mass changes in the HMA area. 
Measurements are made on the glacier surface (Cogley et al., 2011) and the measured parameters 
usually include density in snow pits, changes in relative surface elevation with a network of stakes, 
and ice velocity and thickness near the glacier front and of frontal position (Gardner et al., 2013). 
However, such detailed measurements can be implemented only on several glaciers of the HMA 
area since the glaciers are far inlands with harsh weather conditions and challenging geography 
(Yao et al., 2012). 

To overcome this challenge, assuming that the snow and ice density in the entire HMA are 
accurately determined, satellite DEM (digital elevation model)-based methods are employed to 
derive glacier mass changes according to repeatedly measured surface elevations from satellite 
altimetry such as ICESat (cloud, and land elevation satellite; Zwally et al., 2002) and CryoSat-2 
(Wang et al., 2015), optical stereo pairs measurement such as ASTER (advanced spaceborne 
thermal emission and reflection radiometer; Brun et al., 2017), and bistatic SAR such as SRTM 
(Rodriguez et al., 2006). However, there are obvious uncertainties for mass changes derived from 
these data (Kááb et al., 2012, 2015; Gardner et al., 2013). Especially for satellite altimetry (Kááb 
et al., 2012; Gardner et al., 2013), the following technique and environmental characteristics can 
impact the measurements and results: (1) the ground tracks cover only about the half of the glaciers; 
(2) the repeated ground tracks are commonly off-set by 1-2 km; (3) there is a large signal footprint 
of about 70 m; (4) the ground tracks are repeated only 2—3 times per year; (5) there can be steep 
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glacier slopes; and (6) there can be extensive debris cover and crevasse fields on top of the ice. 

Another option is the usage of space gravimetry. Time-variable (usually monthly) gravity 
changes were derived by continuous measurements of the Gravity Recovery and Climate 
Experiment (GRACE) twin-satellite mission from 2002 to 2017, which can be used to infer the 
monthly glacier mass changes in the HMA area (Bettadpur, 2012). A GRACE data analysis has 
several advantages compared to the glaciological methods and satellite altimetry mentioned above, 
including (1) a full and homogeneous coverage over the glaciers; (2) monthly resolved gravity 
changes and glacier mass changes; and (3) a long observation time span from March 2002 to June 
2017. In turn, adequate post-processing, mainly filtering, of GRACE data must be performed to 
extract the signal of interest as accurately as possible. Most importantly, NS-oriented stripes, which 
are sub-Nyquist artifacts arising from the oversampling of the Earth's low-frequency geoid in EW 
(east to west) direction (Peidou and Pagiatakis, 2020), must be removed with so-called destriping 
filters. 

Several studies used GRACE gravity data to derive glacier mass changes in the HMA area over 
different time spans (Matsuo and Heki, 2010; Gardner et al., 2013). There are, however, obvious 
differences in the estimates of different studies during approximately the same observation periods. 
For example, Jacob et al. (2012) provided trend rates of glacier mass changes for 2003-2010, which 
are —5.0 (+6.0) Gt/a in Himalayas and Nyaingentanglha, —5.0 (+6.0) Gt/a in Tianshan, and —11.0 
(+9.8) Gt/a in the total HMA glacier area (excluding central glacier areas of the Qinghai-Tibet 
Plateau (QTP), and Qilian Mountains (Qilian)). Despite the slightly different observation period 
from 2003 to 2012, Yi and Sun (2014) showed very different results, which are —20.5 (+3.0) Gt/a, 
—8.4 (+2.4) and —35.0 (+4.6) Gt/a, respectively. They used the destriping filter P4M6 (Chambers, 
2006), which apparently could not effectively reduce the NS-orientedstripe noise in the Stokes 
coefficients of the GRACE gravity data, as found in Figure 2c of Xiang et al. (2017). This 
emphasizes the importance of using an adequate post-processing of GRACE data. 

After the end of the GRACE mission in June 2017, the gravity data solutions have been updated 
and a new, so-called release-6 version was released by different institutions, such as the Institute of 
Geodesy at Graz University of Technology (ITSG), Austria (Mayer-Giirr et al., 2016; Kvas et al., 
2019), the Center of Space Research, USA (Bettadpur, 2012), the German Research Center for 
Geosciences (Dahle, et al., 2012), the Jet Propulsion Laboratory, USA (Watkins and Yuan, 2012), 
the Tongji University (Chen et al., 2015) and Huazhong University of Science and Technology 
(Zhou et al., 2016). Following Wahr et al. (1998) and Xiang et al. (2017), we calculate trend rates 
of gravity changes (in EWH (equivalent water height)) from these different new gravity solutions 
in the HMA and adjacent areas, and find that the differences among the results for the models are 
far less than those from different destriping filters (not shown). Therefore, it seems reasonable to 
estimate the glacier mass changes in the HMA area by employing only one of these new GRACE 
gravity solutions but with an effective destriping filter. 

In this study, we investigate the glacier mass changes in the HMA area derived from the ITSG 
release-6 solution together with some meteorological reanalysis data and analyze differences 
between our results and those from previous studies. Due to the lower quality of the Stokes 
coefficients from the solutions near the end of the GRACE mission, the observation time span is 
tailored to April 2002 to August 2016. 


2 Data and methods 


2.1 GRACE gravity data 


The ITSG release-6 solution (i.e., ITSG-Grace2018; Kvas et al., 2019) is available as monthly 
Stokes coefficients from degree 2 to 96 from April 2002 to August 2016 (http://icgem.gfz- 
potsdam.de/series). Due to large errors for the higher d/o Stokes coefficients, the maximum d/o for 
the Stokes coefficients is set to 90. The degree 1 Stokes coefficients from Swenson's GRACE-OBP 
approach (Swenson et al., 2008) are available from ftp://podaac.jpl.nasa.gov/allData/tellus/L2/ 
degree_1/deg1_coef.txt. The C20 Stokes coefficients with large uncertainties (Chen et al., 2016) 
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should be replaced with estimates from satellite laser ranging (Cheng and Tapley, 2004; Cheng et 
al., 2013). The data are from https://podaac-tools.jpl.nasa.gov/drive/files/allData/grace/docs/TN- 
11_C20_SLR.txt. The glacial isostatic adjustment effects are corrected from the monthly Stokes 
coefficients of the ITSG model with the GIA model ICE-6G_C(VM5a) (Peltier et al., 2015). 

According to Xiang et al. (2017), the S&WP2M8 destriping filter (Swenson and Wahr, 2006) is 
recommended since trend rates of mass changes in the HMA and adjacent areas derived from 
different GRACE release-5 solutions agree to the increasing lake level and declining glacier mass 
trends as observed by ICESat data (Gardner et al., 2013; Zhang et al., 2013) when this filter is 
applied. For the S&WP2M68 filter, a degree 2 polynomial is used in a window dependent upon the 
order, fitted to all the Stokes coefficients of the same parity in degrees >8 for a certain order (28), 
and the coefficients used in fitting are subtracted by the fitting values so that the striping errors can 
be reduced in the coefficients within the window. 

Moreover, the effects of soil moisture changes have to be further removed from the monthly 
Stokes coefficients. For this purpose, average grid soil moisture data from two land surface models, 
namely Noah and Variable Infiltration Capacity of the Global Land Data Assimilation System 
(Rodell et al., 2004), are decomposed into spherical harmonic (SH) series using Equations 3-5 of 
Xiang et al. (2016). The calculated SH coefficients can be transformed into the related Stokes 
coefficients according to Equations 1 and 2 of Xiang et al. (2016). Then, these Stokes coefficients 
can be removed from the monthly destriped Stokes coefficients of the ITSG solution. 


2.2 Computation of mass changes 


Firstly, we divide the HMA and adjacent areas into a series of mascon regions with a total number 
of 700, which include 14 irregular glacier mascons based on the Randolph Glacier Inventory 6.0 
(RGI Consortium, 2017), regular 2°x2° mascons in the other regions and unregular mascons in 
their buffer zones. Note that the 14 irregular glacier mascons are allowed to be slightly larger than 
specified in the inventory since RGI 6.0 did not consider both the retreat of glaciers and glacier 
lake expansion from 2001 to 2016. 

Secondly, we remove the average during the observation time span from the Stokes coefficients 
from the ITSG solution and denote the results as dcim and Osim. Supposing that each mascon has a 
uniformly distributed unit EWH, which can be decomposed into SH series, and that the SH 


coefficients Aci, and AS, (i=1, 2, ... , M) can be calculated with Equations 3-5 of Xiang et al. 


(2016), then according to Jacob et al. (2012), the realistic EWH (AEWH)) for the 700 mascons (i=1, 
2, ..., M) can be expressed by, 


i Pave OM paima 21 +1 1 2 i i 
AEWH' = a Za a en Y. Wi (ÓC y ACim + ÓS y AS y)» (1) 
Where a is the Earth's radius; k; is the degree | elastic load Love number for potential perturbation 
(Wang et al., 2012); pave is the average density of the earth (5.517 g/cm); pw is the water density 
(1 g/cm?); wi is the degree coefficient of Jekeli's Gaussian averaging function for Legendre 


expansion (Jekeli, 1981); A is the element of the inverse of matrix 4 while elements of A are 
given by 
Imax 1 2 ; i ; i 
Aj S ae > W; (Ach ACI + ASh AS m) , (2) 
given for a Gaussian filter with an averaging radius of 220 km, the summation of Equations 1 and 
2 is truncated at Imax. 

Thirdly, the monthly mass changes as denoted above by EWH for each mascon are fitted to a 
trend, an annual, a semi-annual and a 161-days aliasing variation (from the S2 semidiurnal solar 
tide) by linear regression, and are corrected with the S2 aliasing signals. 

Furthermore, assuming that the EWH trend estimate at each mascon is disturbed by the 
independent regression error and the noise errors (i.e., measurement errors and remaining stripe 
noises), the uncertainty of the estimated EWH trend can include twice the regression standard 
deviation and the effect of the noise errors. For a target mascon, the effect of the noise errors can 
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be estimated by using the ITSG solution-derived mass change trend magnitude on the ocean at the 
same latitude as the mascon, but excluding ocean area 800 km off the continents to avoid land 
signal leakage (Chen et al., 2009; Longuevergne et al., 2010; Farinotti et al., 2015). 


2.3 Meteorological reanalysis data 


Since some meteorological physical quantities, particularly total net precipitation (net snowfall plus 
net rainfall), net snowfall and total net radiation, could be correlated with the glacier mass changes, 
they can be used for comparison with the estimated glacier mass changes in the selected glacier 
areas. In this study, the monthly total net precipitation is derived from monthly total precipitation 
(snowfall plus rainfall), total evaporation and total runoff data. The monthly net snowfall is 
calculated from monthly snowfall, snow evaporation and snowmelt data. The total net radiation is 
derived from surface net solar radiation and surface net thermal radiation data, respectively. We 
further compare our results to monthly air temperature, which is directly from the monthly 
temperature data of air at 2 m above the surface of land. 

We apply the ERA5-Land reanalysis dataset (Muñoz-Sabater, 2019; Hersbach et al., 2020), 
which combines model data with observations using the laws of physics. The horizontal resolution 
of the dataset is 0.1°x0.1° (native resolution is 9 km). This dataset contains all meteorological data 
used in this study and is available from https://cds.climate.copernicus.eu/cdsapp#!/home. 


3 Results 


3.1 Glacier mass changes 


The derived trend rates of gravity changes expressed by EWH from monthly Stokes coefficients 
(Fig. 1) show that the NS-oriented stripes are effectively suppressed. The negative signals along 
the mountains could be ad hoc related to glacier melting trends. The positive signals in the central 
QTP are likely due to the rising trend of lake water level (e.g., lakes in the southern Qiangtang 
Basin, Zhang et al., 2013) and increasing trend of groundwater in the Three Rivers Source area 
(Xiang et al., 2016). There are also significant negative signals found in Northwest India and the 
Bengal Basin, which are due to the excessive use of groundwater for agricultural irrigation (Rodell 
et al. 2009). 

In order to reduce the signal leakage effects among different glaciers and between glaciers and 
non-glacier regions, monthly mass changes can be simultaneously calculated using Equation 1 at 
the 700 mascons. We focus on the 14 glacier mascons, in which mascon 1 covers the glaciers in 
Nyaingentanglha, mascons 2 to 3 those in Eastern Himalayas, mascons 4 to 5 in Western 
Himalayas, mascon 6 in Karakoram, mascon 7 in Hindu Kush, mascon 8 in Pamirs, mascons 9 to 
12 in Tianshan, mascon 13 in West Kunlun, and mascon 14 in Qilian, respectively. 

We show the trend rates of monthly mass changes at the HMA mascons in Figure 2 and the time 
series for each of the 14 glacier mascons in Figure 3. The values of the mass change trend rates at 
the 14 glacier mascons and some mascon combinations, expressed by both EWH per year and mass 
budget per year (Fig. 3; Table 1). Although results are also available for non-glacier mascons, which 
may show, for example, the declining changes of groundwater storage in Northwest India in Figure 
2, they are not discussed in this study. 

The glaciers experience mass loss except for mascons 6 (Karakoram) and 13 (West Kunlun). As 
can be seen in Table 1, the largest ice melting trends appear at mascon 1 (Nyaingentanglha) and 
mascons 2 to 3 (Eastern Himalayas) with rates of —1.08 (+0.14) m/a (—7.02 (+0.94) Gt/a) and -1.05 
(+0.12) m/a (6.73 (+0.78) Gt/a), respectively. Large ice melting trends are found at mascons 4 to 
5 (Western Himalayas) and mascons 9 to 12 (Tian Shan) with rates of —0.36 (+0.07) m/a (43.46 
(+0.66) Gt/a) and —0.48 (+0.06) m/a) (55.10 (+0.68) Gt/a). Small but well-determined glacier 
melting trends are found at mascon 7 (Hindu Kush) and mascon 14 (Qilian) with rates of -0.23 
(+0.11) m/a (0.97 (+0.47) Gt/a) and —0.65 (+0.28) m/a (20.76 (+0.33) Gt/a). A slight glacier 
melting trend of —0.01 (+0.08) m/a (-0.08 (+0.83) Gt/a) results at mascon 8 (Pamirs), in which the 
estimated signal value is far less than the uncertainty. In turn, weak mass gain trends of 0.06 (+0.03) 
m/a (1.19 (+0.55) Gt/a) and 0.11 (+0.05) m/a (0.77 (+0.37) Gt/a) are found at mascon 6 


202104.00095v1 


chinaXiv 


ChinaXivA ERAT! 


XIANG Longwei et al.: Glacier mass balance in High Mountain Asia inferred from a GRACE... 


(Karakoram) and mascon 13 (West Kunlun), respectively. In total, the HMA glaciers are found to 
have a mass loss trend of —22.17 (+1.96) Gt/a. 

As indicated in the mass trend results in Table 1, the glacier areas in the Karakoram, West Kunlun 
and Pamir in the northwest of HMA form a so-called 'Karakoram anomaly’ indeed shows an almost 
zero or positive mass balance (de Kok et al., 2018; Farinotti et al., 2020). The other areas usually 
show a significant mass loss trend (Yao et al., 2012). 
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Fig. 2 Trend rates of the monthly mass changes (in EWH) at the mascons in the HMA area and adjacent regions 
during the observation time for the period April 2002 to August 2016, inverted from the GRACE release-6 ITSG 
solution. In the inversion process, the Stokes coefficients are destriped with the S&W P2M8 filter and smoothed 
by a Gaussian filter with averaging radius of 220 km. Monthly mass changes at the 14 glacier mascons are shown 
in Figure 3. 


Table 1 Trend rates of glacier mass changes at the 14 glacier mascons and selected mascon combinations 


Mascon Mascon area Trend (Gt/a) Area (km?) Trend (m/a w. e.) 

1 Nyaingentanglha Mountains —7.02+0.94 6498.14 -1.08+0.14 

2 -2.72+0.56 1909.14 -1.42+0.29 
Eastern Himalayas 

3 —4.01+0.55 4527.44 —0.89+0.12 

4 —1.10+0.44 4128.77 —0.27+0.11 
Western Himalayas 

5 -2.36+0.50 5413.30 —0.44+0.09 

6 Karakoram Mountains 1.19+0.55 20638.73 0.06+0.03 

7 Hindu Kush Mountains —0.97+0.47 4299.61 -0.23+0.11 

8 Pamirs —0.08+0.83 10580.46 —0.01+0.08 

9 -1.06+0.34 7106.08 —0.15+0.05 

10 —0.89+0.43 952.15 —0.93+0.46 
Tianshan Mountains 

11 —2.74+0.38 2166.41 —1.26+0.17 

12 —0.41+0.13 415.29 —0.99+0.32 

13 Western Kunlun Mountains 0.7740.37 7159.88 0.11+0.05 

14 Qilian Mountains -0.76+0.33 1176.82 —0.65+0.28 

Total Total HMA glaciers -22.17+1.96 76972.22 —0.29+0.03 

2+3 Eastern Himalayas —6.73+0.78 6436.58 —1.05+0.12 

4+5 Western Himalayas —3.46+0.66 9542.07 —0.36+0.07 

9+10+11+12 Tianshan Mountains —5.10+0.68 10639.93 —0.48+0.06 


Note: HMA, High Mountain Asia. 
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Fig. 3 GRACE-derived monthly mass changes at the 14 glacier mascons 
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4 Discussion 


4.1 Comparison with previous estimates 


In Tables 2 and 3, the trend rates of the GRACE-derived monthly glacier mass changes are 
calculated for different glacier mascons and glacier mascon combinations during different 
observation time spans, in order to compare them with those of previous studies having the same 
observation time spans. Our results can match those from ICESat DEM and ASTER DEM analysis 
at most glacier mascons or glacier mascon combinations. However, some results show larger 
differences (Table 2). 


Table 2 Comparisons of trend rates of monthly glacier mass changes between this study and previous studies at 
the ten selected glacier mascons during the same observation time span 


Trend rates of monthly glacier mass change (Gt/a) 


Mascon ICESat ASTER GRACE 
2003-2008 2003-2009 2003-2009 2000-2016 2003-2008 2003-2009 2002-2016 
1 —8.0+1.7 —4,543.3 -8.3+4.33 —4.0+1.5 —9.12+2.70 —8.67+2.29 —7.02+0.94 
2 -3.1+0.6 —1.0+0.5 —0.82+1.66 —2.72+0.56 
—5.5+1.6 —4.4+2.1 —6.18+2.01 
3 —3.2+0.9 —1.6+1.0 —4.12+1.72 -4.01+0.55 
4 —3.2+0.7 —3.1+1.8 —0.9+0.6 —1.6+0.4 —3.14+1.28 —2,.38+1.09 —1.10+0.44 
5 —4.7+1.1 —4.4+1.6 - —2.9+0.7 —3.01+1.47 —3.19+1.24 —2.36+0.50 
6 —2.1+1.3 - -0.5+1.2 1.75+1.66 1.19+0.55 
—2.6+4.4 -1.00+1.86 
7 —2.7+0.6 - —0.6+0.4 -4.36+1.43 —0.97+0.47 
8 —3.1+0.9 -2.1+4.1 - —0.7+0.5 —2.50+2.05 —0.72+1.89 —0.08+0.83 
13 0.6+0.9 1.5+1.7 0.2+1.6 1.4+0.8 1.21+1.09 0.28+0.95 0.77+0.37 
14 - —0.6+0.8 - - - —0.68+0.83 - 
Source Kääb et al. Gardner etal. Neckel et al. Brun et al. This stud 
(2015) (2013) (2014) (2017) y 


Note: Larger differences in previous results are underlined. -, no result available. 
Table 3 Comparisons of trend rates of monthly glacier mass changes between this study and previous studies for 
selected glacier mascon combinations 
Trend rates of monthly glacier mass change (Gt/a) 
Mascon ICESat GRACE ASTER GRACE 
2003-2009 2003-2010 2003-2012 2000-2016 2003-2009 2003-2010 2003-2012 2002-2016 
1+2+3+4+5 -17.5+4.4 —5.0+6.0 —20.543.0  -11.1+2.0 -20.42+3.46 -21.88+2.81 -19.00+2.64 -17.21+1.39 


a7 3.2462 1050 #25 416 -1.444281 1.914213 1974179 9914116 
9+10+11+12 -7.5434 5.0460 8442.4 -3.042.2 -7.7141.75 -5.41+1.33 -4.38+1.18 -5.10+0.68 
Total -28.2+8.3 -11.049.8  -35.0+4.6 -14.543.4 -29.57+4.79 -25.48+3.77 -21.41+3.40 -21.40+1.93 
Source Gardner et al. Jacob etal. Yiand Sun Brun etal. This stud 
(2013) (2012) (2014) (2017) y 


As can be seen in Table 2, Kááb et al. (2015) provide a value of —3.1 (+0.6) Gt/a for mascon 2 
(eastern half of the Eastern Himalayas) from 2003 to 2008 while our result of -0.82 (+1.66) Gt/a is 
much smaller, i.e., only about one fourth. The combined value for mascons 2 and 3 together (-6.3 
(+1.1) Gt/a) agrees to ours (-4.94 (+2.39) Gt/a) within the errors as our estimate for mascon 3 is 
larger than the one from Kááb et al. (2015). This may point to improvements in the leakage 
correction. 

Gardner et al. (2013) listed the value of —4.5 (43.3) Gt/a for mascon 1 (Nyaingentanglha) from 
2003 to 2009 while our result is -8.67 (+2.29) Gt/a, almost twice as much. As our value for mascon 
1 agrees very well to those of Kááb et al. (2015) and Neckel et al. (2014), the value of Gardner et 
al. (2013) appears to be an outlier which is suggested to be reinvestigated. Neckel et al. (2014) 
calculated —4.4 (+2.1) Gt/a for mascons 2 to 3 (Eastern Himalayas) from 2003 to 2009 while our 
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result of —6.18 (+2.01) Gt/a is 40% larger for the same time span. However, despite these 
differences these values still agree within their errors. In particular, Gardner et al. (2013) generally 
provide larger errors than the other studies or our GRACE-based results. 

Brun et al. (2017) determined —4.0 (+1.5) Gt/a for mascon 1 (Nyaingentanglha) from 2000 to 
2016 while our result of —7.02 (+0.94) Gt/a from 2002 to 2016 is 75% larger. Except for mascons 
1-3, there is good agreement between our results and the ones from Brun et al. (2017). For mascons 
2 and 3, we speculate that an improved leakage correction in the GRACE data could reduce the 
differences, but it should also be investigated if the ASTER DEM data are satisfactory. Similarly, 
for mascon 1 the apparent difference might relate to limitations in the ASTER DEM data, 
especially, because our values for the shorter time spans agree very well to two other studies. 

We also evaluate how our results perform in comparison to combined regional results (Table 3). 
Gardner et al. (2013) listed a value of —17.5 (+4.4) Gt/a for mascons 1 to 5 (Nyaingentanglha and 
Himalayas) from 2003 to 2009 while our result of —20.42 (+3.46) Gt/a is slightly larger. Brun et al. 
(2017) gave values of —11.10 (+2.0) Gt/a for mascons 1 to 5 (Nyaingentanglha and Himalayas) and 
-3.0 (+2.2) Gt/a for mascons 9 to 12 (Tianshan) from 2000 to 2016 while our results, -17.21 (+1.39) 
Gt/a and -5.10 (+0.68) Gt/a, respectively, but from 2002 to 2016, are both larger. We think that the 
partly significant differences to those previous studies based on the ICESat and ASTER DEMs 
analysis are due to limitations such as coverage, spatial-temporal resolution and noise level of the 
other satellite data (Gardner et al., 2013; Brun et al., 2017). According to Farinotti et al. (2015), the 
trend rates of mass changes for the glacier mascons 9 to 12 from 2003 to 2009 were estimated to 
be 6.2 (+2.8) Gt/a by glaciological modelling, which is close to our estimate (-5.41 (+1.33) Gt/a). 
Based on Gardner et al. (2013), all the glaciers together had a mass loss trend of —28.2 (+8.3) Gt/a 
from 2003 to 2009, which agrees very well to our estimate of -29.57 (+4.79) Gt/a. Note that ASTER 
DEMs (digital elevation models) compute the mass balance of 92% of the glaciated area while the 
remaining 8% including small glaciers were not processed (Brun et al., 2017), and thus might miss 
parts of the mass changes. 

Table 3 also shows the comparison of our GRACE-based results to two previous GRACE-based 
studies (Jacob et al., 2012; Yi and Sun, 2014). Large differences are also found here. For example, 
Jacob et al. (2012) determined -5.0 (+6.0) Gt/a for mascons 1 to 5 (Nyaingentanglha and 
Himalayas) from 2003 to 2010 while our result of —21.88 (+2.81) Gt/a is more than 4 times larger, 
but with much smaller uncertainty. Yi & Sun (2014) gave the values of —6.1 (+2.5) Gt/a for mascons 
6 to 8 (Karakoram, Hindu Kush and Pamirs) and —8.4 (+2.4) Gt/a for mascons 9 to 12 (Tianshan) 
from 2003 to 2012 while our results are with 1.97 (+1.79) Gt/a and -4.38 (+1.18) Gt/a, respectively, 
much smaller. These large differences between our results and those from two previous GRACE 
based studies are referred to different methods to suppress the NS-oriented stripe noise and 
measurement errors, as well as different GRACE solutions. 


4.2 Comparison with changes in meteorological quantities 


In Figure 3, the time series of monthly glacier mass changes show manifold time-variable features. 
We therefore calculate interannual (Figs. 4 and 5) and average seasonal mass changes (Fig. 6) and 
compare them with corresponding total precipitation, snowfall, surface air temperature and total 
radiation data at nine selected glacier areas approximately covering all the HMA glaciers. Note that 
total precipitation contains snowfall and rainfall. The latter is detrimental to accumulation as it 
carries heat. 


4.2.1 Interannual changes 


An increase/decrease in snowfall could induce more/less glacier mass accumulation. Also, an 
increase/decrease in summer air temperature could possibly cause more/less glacier mass ablation. 
The changes of these two quantities correspond to partial glacier mass changes or mass balances 
(Yao et al., 2012; Farinotti et al., 2020), thus they are tentatively used here to validate the observed 
glacier mass changes. Additionally, we compare total precipitation measurements to our results. Our 
GRACE-derived monthly glacier mass changes, and ERA5-Land monthly total precipitation and 
snowfall are averaged over 12 months to get annual averaged results for each quantity in nine selected 
glacier mascon areas. The ERA5-Land summer temperatures are averaged over the 3 summer months 
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of a year (June, July and August) to get three-summer-month averaged temperatures. These annual 
results are shown for comparison in Figure 4. However, there are no significant trends found for the 
total precipitation, snowfall and the difference between total precipitation and snowfall, i.e., rainfall. 
The maximal magnitudes for the total precipitation and snowfall are not larger than 0.08 and 0.07 
Gt/a, respectively, which are far less than the observed glacier mass changes (Table 1). In addition, 
the increasing temperature trends could help explain the mass loss trends for mascon 1 
(Nyaingentanglha), mascons 2 to 3 (Eastern Himalayas), mascon 7 (Hindu Kush), mascon 8 (Pamirs), 
mascons 9 to 12 (Tianshan) and mascon 14 (Qilian), but fail to do so for mascon 13 (Western Kunlun) 
that shows a mass gain trend. A decreasing temperature trend (—0.015°C/a) could correspond to a 
mass gain trend for mascon 6 (Karakoram), but it (-0.043°C/a) cannot support a mass loss trend in 
mascons 4 to 5 (Western Himalayas). 
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Fig. 4 Annual averaged glacier mass changes showing interannual glacier mass changes in the nine selected 
glacier areas compared with the annual changes in total precipitation and snowfall, and three-summer-month 
averaged surface air temperatures, respectively. The trend rates of air temperature are also shown for the nine glacier 
areas. 


In view of the unsatisfactory comparisons above, we additionally compare our annual glacier 
mass changes, as suggested in recent studies (de Kok et al., 2018; Farinotti et al., 2020), with 
changes in total net precipitation and net snowfall, and total net radiation energy (Fig. 5). The net 
snowfall in the HMA is directly related to the glacier mass change (Brock et al., 2010). Not that the 
discrepancy between total net precipitation and net snowfall is actually net rainfall, which may 
partly directly contribute to the glacier mass change. And changes in net rainfall could also imply 
the heat effect on the glacier mass change. Since the total energy balance largely controls glacier 
ablation, changes in total net radiation directly impact the glacier mass balance (Brock et al., 2010). 
The results of three meteorological quantities are accumulated in order to make a precise 
comparison with the already accumulated mass changes derived from GRACE data. Note that 
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before accumulating, the 14-a average is removed from the monthly net radiation energy. Still, the 
mass change results cannot be fully explained with these data. First, total net precipitation and net 
snowfall usually have the same trend signs except for mascons 9 to 12 (Tianhan) and mascon 14 
(Qilian). Both could thus explain the mass gain trends at mascon 6 (Karakoram), mascon 8 (Pamir) 
and mascon 13 (Western Kunlun). The decreasing trend in net snowfall together with the increasing 
trend in net rainfall causing heat effect could be related to the mass loss trend for mascons 9 to 12 
(Tianshan) and mascon 14 (Qilian). Total net radiation energy shows an increasing trend for all 
glacier areas, which could imply a mass loss trend for mascon 1 (Nyaingentanglha), mascons 2 to 
3 (Eastern Himalayas), mascons 4 to 5 (Western Himalayas), mascon 7 (Hindu Kush), mascons 9 
to 12 (Tianshan) and mascon 14 (Qilian). It is confirmed that an increase/decrease in non-trend 
total net radiation energy can usually explain a decrease/increase in glacier non-trend mass change 
for each glacier area. In addition, because of the heat effect, the increasing trend in net rainfall 
could also cause a mass loss trend for mascon 1 (Nyaingentanglha), and mascons 2 to 3 (Eastern 
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Fig. 5 Annual averaged glacier mass changes in the nine glacier areas are compared with the annual changes in 
accumulative total net precipitation, accumulative net snowfall and accumulative total net radiation energy 


For the glacier areas of the so-called 'Karakoram anomaly’, the mass gain trends are usually 
thought to be mainly due to an increase in the strength and frequency of the westerlies-dominated 
total precipitation or snowfall during the last several decades (Cannon et al., 2015), and they are 
also supported by an increase in summer snowfall caused by the increased land irrigation in the 
lowlands (de Kok et al., 2018). Other mechanisms are not introduced (de Kok et al., 2018; Farinotti 
et al., 2020). For other glacier areas, the mass loss trends could be mainly due to the weakening 
Indian monsoon bringing reduced total precipitation. However, as stated above, there are no 
obvious trends in total precipitation and snowfall during the 14-a observation time that could 
support this view. Although we also tested accumulative total net precipitation, net snowfall data 
and net radiation energy data, they could only explain mass changes for parts of the glacier regions. 
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A possible reason for this is that the used meteorological reanalyses data have larger uncertainties, 
especially for the trend estimates (Palazzi,et al., 2013). 


2.40, 0.60 20 
1.20 0.30 10 
0.00 0.00 0 
-1.20 0.304 -10 
2.40 0.60! -20 600 
Feb Apr Jun Aug Oct Dec Feb Apr Jun Aug Oct Dec 
4.00 30 600 4.00, 30 600 
Mascons 4-5 nn > 
2.00 15 A 300 5 2.0013 
0.00 0 \ o $ 0.00} A 
-2.00} _ -1 -300 -2.00 £ 
5 4.001 9-3 -600 -E-4.00l 2 ~ 
= Feb A J A Oct D a 3 Feb A; J A Oct D 
= 0.80 8, 30 e pr Jun Aug Oct Dec 605 É 2.40, E 20 el pr Jun Aug Oct Dec soo 
3 0.40} £ Fa 
A S 3 
‘5 0.00 A LO 
2 £ 3 
£ -0.40} =-1 3 r 
Ž -0.80L 3-30 600 
5 2 Feb Apr Jun Aug Oct Dec ea Feb Apr Jun Aug Oct Dec 
1.40, Y 15 600 0.20, 600 
1 Mascon 13 
0.70 oio 5 SNL 300 
0.00 0.00: 0 ; 0 
-0.70 -0.10 -5 Vv -300 
1.40 0.20. -10 n s a A n +600 
Feb Apr Jun Aug Oct Dec 
0.06 Month 
0.03 — Glacier mass change 
0.00 — Accumulated net precipitation change 
—0.03 — Accumulated net snowfall change 
—0.06 ~ Accumulated net radiation change 


Feb Apr Jun Aug Oct Dec 
Month 


Fig.6 14-a averaged glacier mass changes for the 12 individual months of one year showing an average seasonal 
mass change in the nine glacier areas, in comparison to corresponding accumulative total net precipitations, 
accumulative net snowfalls and accumulative total net radiation energies, respectively. 


4.2.2 Average seasonal changes 


We further calculate a 14-a average for each of the 12 months of one year from monthly glacier 
mass changes, ERA5-Land monthly accumulative total net precipitations, accumulative net 
snowfalls and accumulative total radiation energy results, respectively. Before accumulating, the 
14-a average values are removed from the monthly net radiation energy results. Moreover, before 
averaging, we remove the corresponding interannual signals and trend signals for the four 
quantities. The comparison results are shown in Figure 6. Seasonal mass changes peak in May- 
July, and show a valley in October-November for mascon 1 (Nyaingentanglha), mascons 4 to 5 
(Western Himalayas), mascon 6 (Karakoram) and mascon 13 (Western Kunlun), respectively. 
Similar seasonal mass changes are also found for mascon 7 (Hindu Kush) and mascon 8 (Pamir), 
but with a peak in April, and a valley in August-September. These seasonal mass changes except 
for mascon 13 are found to agree with the changes in accumulative total net precipitation and net 
snowfall and are anti-correlated to accumulative total radiation energy. The seasonal changes in 
accumulative net snowfall almost account for all those in the total net precipitation. For mascon 13 
(Western Kunlun), the seasonal mass change could be explained by both increased accumulative 
net snowfall in spring and increased accumulative net rainfall in summer. However, for other glacier 
areas, i.e., mascons 2 to 3 (Eastern Himalayas), mascons 9 to 12 (Tianshan) and mascon 14 (Qilian), 
there are more peaks and valleys for the seasonal mass changes resulting a more complicated 
pattern, which can obviously not be explained by the changes in the three meteorological quantities. 
The reason for this is unclear and deserves further analysis. 
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5 Conclusions 


We provided estimates of glacier mass changes in the HMA area from April 2002 to August 2016 
by employing a new GRACE solution. The total mass loss trend of -22.17 (+1.96) Gt/a is found 
for the HMA glaciers. The largest glacier mass losses yield in the Nyaingentanglha and the Eastern 
Himalayas with rates of —7.02 (+0.94) and -6.73 (+0.78) Gt/a, respectively. The glacier areas within 
the so-called 'Karakoram anomaly' have a nearly zero mass balance in the Pamirs and a slight mass 
gain trend in the Karakoram and West Kunlun with rates of 1.19 (+0.55) and 0.77 (+0.37) Gt/a, 
respectively. 

Our estimates of glacier mass change trends can match the results from previous studies 
according to the ICESat and ASTER DEM analysis in most of the selected glacier areas where 
results are available over the same time periods for comparison. Larger differences in some glacier 
areas could be due to their limitations on the coverage, spatial-temporal resolution and noise level 
of data. However, our results are found to be very different from those of two previous GRACE- 
based studies, possibly due to our different post-processing techniques and the new GRACE 
solution used by us. 

The 14-a estimates of glacier mass change for the nine selected glacier areas explicitly show 
different interannual mass changes and average seasonal mass changes, which can be partly 
supported by total net precipitation and net snowfall data, but mostly seem to correlate with total 
net radiation energy data. Almost all non-trend interannual mass changes and most seasonal mass 
changes can be explained by total net radiation energy data. The mass loss trends could be partly 
due to heat effect from increased net rainfall in some glacier areas such as in Tianshan, Qilian, 
Nyaingentanglha and Eastern Himalayas. 

The new glacier mass change results in this study based on the latest GRACE solution can help 
to serve the utilization of water resources in the HMA and adjacent areas and contribute to global 
change studies. 
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